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Scaling of the critical slip distance in granular layers 

Takahiro Hatano 1 



We investigate the nature of friction in granular layers by 
means of numerical simulation focusing on the critical slip 
distance, over which the system relaxes to a new stationary 
state. Analyzing a transient process in which the sliding 
velocity is instantaneously changed, we find that the criti- 
cal slip distance is proportional to the sliding velocity. We 
thus define the relaxation time, which is independent of the 
sliding velocity. It is found that the relaxation time is pro- 
portional to the layer thickness and inversely proportional 
to the square root of the pressure. An evolution law for 
the relaxation process is proposed, which does not contain 
any length constants describing the surface geometry but 
the relaxation time of the bulk granular matter. As a re- 
sult, the critical slip distance is scaled with a typical length 
scale of a system. It is proportional to the layer thickness in 
an instantaneous velocity change experiment, whereas it is 
scaled with the total slip distance in a spring-block system 
on granular layers. 



1. Introduction 

A natural fault has the cataclasite core zone, along which 
shear deformation concentrates (e.g. Engelder [1974]; Scholz 
[1987]). Rheology of these granular matters thus provides 
us an important insight in considering the nature of friction 
on faults from a microscopic point of view. Unfortunately, 
to this date, our understanding of the rheological proper- 
ties of granular matter is still poor except for dilute flow 
to which the kinetic theory of gases can apply (e.g. Garzo 
and Dufty [1999] and references therein.) Thus, a computa- 
tional approach has played a considerable role in investigat- 
ing dense granular rheology to propose some constitutive 
laws for stationary shear flow (e.g. Aharonov and Sparks 
[2002]; Hazzard and Mair [2003]; GDR MiDi [2004]; da Cruz 
et al. [2005]; J op et al. [2004]; Hatano [2007]). However, a 
nonstationary state is still a frontier in the sense that we do 
not have any constitutive laws for transient processes. 

The description of transient states is particularly impor- 
tant in the context of seismology because an earthquake is 
essentially a nonstationary process. An important quantity 
is the critical slip distance, over which a fault looses its fric- 
tional strength with the coseismic slip [Ida, 1973], because it 
determines the maximum acceleration of the seismic ground 
motion [Aki, 1987] as well as the rupture nucleation process 
(e.g. Ohnaka [2000]). However, regardless of its importance, 
we still can not explain the critical slip distance ranging from 
10" 1 to 1 m, which is obtained by the seismic inversion [Ide 
and Takeo, 1997]. It is rather paradoxical that the critical 
slip distance obtained in a typical experiment is of the order 
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of 10~ 5 m [Dieterich, 1979; Scholz, 2002]. Understanding 
the physics that determines the critical slip distance to ex- 
plain the wide gap between a natural fault and a laboratory 
is thus a central problem in seismology [Marone and Kilgore, 
1998; Ohnaka and Shen, 1999; Ohnaka, 2003]. 

In this letter, by means of numerical simulation on 
sheared granular layers, we obtain a constitutive law that 
describes a nonstationary process as well as a stationary 
state. Using this constitutive law together with dimensional 
analysis, we propose a new interpretation on the scale de- 
pendence of the critical slip distance. 

2. Model 

In the following we describe the computational model of 
granular layers. Each grain is assumed to be sphere. The 
interaction between grains is described by the discrete el- 
ement method (DEM) [Cundall and Strack , 1979], which 
is the standard model used in powder engineering and soil 
mechanics. Consider a grain i of radius Ri located at n 
with the translational velocity Vi and the angular velocity 
Cli. This grain interacts with another grain j when they are 
in contact; i.e. |ry| < Ri + Rj, where ry = r; — r^. The 
interaction consists of two kinds of forces, which are nor- 
mal and transverse to ry , respectively. Introducing the unit 
normal vector ny = r y/l r ij'l> the normal force acting on i, 
which is denoted by F'j', is given by [khij + £ny ■ ry] ny, 
where hij = Ri + Rj — |ry|. In order to define the trans- 
verse force, we utilize the relative tangential velocity 
defined by ry — n t j • ry + (Rifli + RjClj)/(Ri + Rj) x ry 
and introduce the relative tangential displacement vector 
AVV = J dtv^ . The subscript in the integral indicates 
that the integral is performed only when the contact is 
rolling; i.e., fct|Ay-| < PelF^'l or Ay • «y < 0. Other- 
wise, the contact is said to be sliding. The expression of 
the tangential force depends on the state of the contact: 
l«e|Fy |Vj*Vl v y I f° r sliding contact and feA^' for rolling 
contact. 

We consider a bidisperse system in order to avoid crys- 
tallization, where the diameters of the constituent particles 
are 0.7a! and 1.0(2, respectively. The number of each grain is 
the same. For simplicity, we assume that the mass of these 
grains is the same, which is denoted by m. The dimensions 
of the system are LxLxH, where we use periodic boundary 
conditions along the x and the y axes. In the z direction, 
there exist two rigid walls that consist of the larger particles. 
One of the walls is displaced along the x axis at constant 
velocity V to realize plain shear flow, where the velocity gra- 
dient is formed in the z direction. These grains interact with 
the bulk grains via the force described above. This wall is 
also allowed to move along the z axis so that the pressure 
is kept constant at P, while it is immobile along the y axis. 
Namely, the z dimension of the system, denoted by H, is a 
dynamic variable. The equation of motion of the wall along 
the z axis is given as MH — F z — PL 2 , where M denotes 
the mass of the wall and F z is the repulsive force given by 
the grains. We check that the mass of the wall does not 
influence the result. Hereafter we set M = 100m. We also 
confirm that the x and y dimensions do not influence the 
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result by comparing two systems: L — Wd and 20d. We 
thus adopt L = 10d to save the computational time. The 
traction acting on the moving wall is monitored, which is 
denoted by F x , so that the friction coefficient of the system 
is defined as = F x /L? /P. We do not consider gravity 
here. 

For de-dimensionalization, we set d = 1, k = 1, and 
m — 1. The other parameters are chosen as kt = k/5, 
[i e — 0.6, and £ = 1. The coefficient of restitution vanishes 
with these parameters, which may be justified in model- 
ing rock powder. However, the general correspondence of 
these parameter values in DEM to the material constants 
in real life is not rigorous. Here we estimate the material 
constants based on the sound velocity and the Y oung's mod- 
ulus, which are roughly estimated as d^Jk/m and k/d, re- 
spectively. (Note that the numerical factors are neglected.) 
Thus, the unit velocity and the unit pressure in DEM are of 
the order of kilometer per second and several tens of Giga- 
pascal, respectively. 

3. Result 

We first prepare a stationary state with the wall velocity 
V\. The stationarity is checked by monitoring the friction 
coefficient, the volume, and the internal velocity gradient. 
Then the wall velocity is instantaneously switched from Vi 
to V 2 . The typical response of the system is shown in Figure 
1, where the friction coefficient relaxes to a new stationary 
value corresponding to the new sliding velocity V 2 . The 
sharp increase of the friction coefficient at the instance of 
the velocity change is due to the steeper velocity gradient 
near the wall. Then this nonlinear velocity gradient relaxes 
to the uniform velocity gradient, which leads to the relax- 
ation of the friction coefficient. The relaxation behaviors of 
the friction coefficient fi and the layer thickness H can be 
fitted by the exponential curve, as is shown in Figures la 
and b. 



H(x) — [l 2 + Qui - 

H(x) = H 2 + (H 1 



/j, 2 )exp(~x/D c ), 
- H 2 )exp(-x/D c 



(1) 

(2) 



where x denotes the slip distance after the velocity switch 
and D c defines the critical slip distance. We confirm that 
the critical slip distance is almost the same for the friction 
coefficient and the layer thickness. We test several values 
of Vi, V 2 , P, and H to find that the choice of Vi does 
not apparently affect the critical slip distance. We thus fix 
V\ = 1 X 10 -5 and change V2, P, and H in the following. 

As is shown in Figure 2, the critical slip distance depends 
on the slip velocity V2 and the normal pressure P. Impor- 
tantly, the critical slip distance is proportional to the slip 
velocity so that the relaxation time can be defined as the 
proportional coefficient. 



D c ~ tV 2 , 



(3) 



where the relaxation time r is independent of the velocity 
(but still depends on the pressure) . In this sense, the relax- 
ation time is more fundamental than the critical slip distance 
in sheared granular layers. This makes a quite contrast to 
conventional experiments on friction of two surfaces, where 
the critical slip distance is independent of the sliding veloc- 
ity (e.g. see Marone [1998]). The discrepancy is due to the 
different physical mechanisms of friction. In the rubbing of 
two surfaces, the area of true contact (asperity) determines 
the friction coefficient so that the critical slip distance is of 
the order of the typical dimension of asperities (e.g. tens 



of micrometers), whereas the internal velocity profile deter- 
mines the friction coefficient in granular layers. 

Then we discuss the nature of the relaxation time. In 
Figure 3a, we find that the relaxation time is inversely pro- 
portional to the square root of the pressure. 



-1/2 



(4) 



Although it has been recognized that the intrinsic time con- 
stant in granular matter should be scaled as equation (4) 
from the viewpoint of dimensional analysis, this relation has 
not been confirmed in a dense system. It is also found that 
the relaxation time is proportional to the layer thickness, as 
is confirmed in Figure 3 a. 



oc H, 



(5) 



which implies that the perturbation propagates into the 
granular layers at the constant velocity. This makes a quite 
contrast to Newtonian fluids, where the velocity field is diffu- 
sive so that the relaxation time is proportional to the square 
of the layer thickness. However, at this point, we cannot de- 
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Figure 1. The relaxation behaviors of the system after 
the velocity change. Here Vi=lx KT 5 , V% = 3 X 10" 3 , 
and P = 1 x 10 -3 . (a) The friction coefficient /i and (b) 
the layer thickness H . The horizontal axes represent the 
slip distance after the velocity change. Symbols denote 
the simulation data, while the solid lines denote the ex- 
ponential curves, Eqs. (1) and (2). (c) Relaxation of the 
internal velocity profile, which is defined as the instanta- 
neous local mean velocity in the x direction. 
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Figure 2. The critical slip distance D c as a function of 
the velocity, V = V 2 . (a) H ~ 6d, (b) H ~ 12d, and (c) 
H ~ 24d. 
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rive Eqs. (4) and (5) from the microscopic principle, i.e., 
the particle dynamics, 

From Eqs. (4) and (5), the relaxation time reads 




where c is a numerical factor. We find c=1.0±0.1in Figure 
3a. Equivalently, from equation (3), the critical slip distance 
can be written as 



where / = V2/H ^Jm/Pd is a nondimensional number 
known as the inertial number (see GDR MiDi [2004] for its 
usefulness in describing granular rheology, while it is origi- 
nally defined by Savage and Hutter [1989].) We can confirm 
the validity of equation (7) in Figure 3b. Note that the crit- 
ical slip distance is proportional to the layer thickness. This 
is consistent with an experiment in which the critical slip 
distance is scaled with the gouge layer thickness [Marone 
and Kilgore, 1998]. 

Then we introduce the evolution law for the transient 
states. The exponential relaxation of the friction coefficient, 
equation (1), implies that the evolution law is a first order 




Figure 3. (a) The pressure dependence of the relaxation 
time, defined by equation (3). Each line denotes r /P 1 * 12 = 
6, 14, 24. (b) The critical slip distance Dc divided by H as 
a function of the inertial number / = Vi /H yj m/Pd mul- 

tiplied by H/d, i.e., V2/ ' yj m/Pd 3 . Thus, this indicates 
that the critical slip distance is proportional to the layer 
thickness. The solid line represents D c /H = V^y/m/Pd 3 . 
The legends are the same as those in Figure 2. 




Figure 4. Unstable sliding of a spring-block system, de- 
scribed by equations (10) and (14), due to the static fric- 
tion, (a) Velocity dependence of the friction coefficient 
during the unstable slip, (b) Temporal behavior of the 
block velocity. 



linear differential equation. Because the critical slip distance 
depends on the sliding velocity, it is more convenient to de- 
scribe the evolution law with respect to time instead of the 
slip distance. The relaxation process is then described by 
the following evolution law. 

A(t) = -r(/i(t) - /*.), (8) 

where /i ss is the stationary friction coefficient, which gener- 
ally depends on the inertial number and other nondimen- 
sional parameters. We do not discuss a stationary constitu- 
tive law here (See, for example, [GDR MiDi, 2004; da Cruz 
et al. , 2005; Jop et al, 2004; Hatano, 2007]). It is essential 
to notice that equation (8) is length-free; i.e., the equation 
does not have any length constants. In the following section, 
we discuss some important consequences of this length-free 
evolution law. 

4. Discussions 

So far we have discussed the nature of the relaxation time 
in a transient process of granular layers and obtained the 
evolution law. In the rest of this paper, we discuss their im- 
portant consequences in application to a more general sit- 
uation. First, using equation (8), we describe the unstable 
sliding of a block on granular layers. The equation of motion 
is given as 

MX = -MQ. 2 (X - X ) - fiN, (9) 

where M is the mass of a block, N is the normal load, MQ 2 
is the stiffness of the spring, and Xo is the equilibrium point 
of the spring. As the temporal evolution of the friction coef- 
ficient /1 is given by equation (8), the block motion is deter- 
mined combining these equations. Before explicitly solving 
Eqs. (8) and (9), we can see the essential property of the 
dynamics by dimensional analysis. Note that equation (9) 
has the time constant Q' 1 and the length constant N /MQ, 2 , 
whereas equation (8) does not have any length constants but 
the time constant r. This means that any lengths defined in 
the resulting dynamics, such as the critical slip distance, is 
scaled by N/MQ. 2 . Namely, the critical slip distance is not 
prescribed by the characteristic length of the microscopic ge- 
ometry but more macroscopic parameters: the normal load 
and the stiffness of the spring. In order to see this more 
explicitly, we solve Eqs. (8) and (9). To this end, it is 
convenient to de-dimensionalize equation (9) as 

£T 2 x — —x — /i, (10) 

where x = (X - X )MO, 2 /N. The initial condition of the 
block motion is given as x(0) = — £t(0) and x(0) = 0. We 
assume that the slip instability is caused by the static fric- 
tion; i.e., fj,(0) — fj, s . This is easily realizable in granular 
layers if the initial state is sufficiently consolidated. Here 
we consider only x(t) > 0, i.e., we solve the block motion 
until it stops again. 

The dynamics of the block is explicitly solvable if we as- 
sume that the dynamic friction coefficient is independent of 
the sliding velocity; i.e., uss = M<2- Then the solution is 
given as 

x{t) = - ^+£2 [C sinQt + costlt + C 2 e- t/T ] - Md,(H) 

fi{t) = (/x s - Md)e"* /T + Hd, (12) 

where C = Qt is the ratio of the two time constants in Eqs. 
(8) and (10). Then the the critical slip distance is apparent 
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from Eqs. (11) and (12). 

N 

D c ~[x(r)+ Ms ] w . (13) 

The numerical factor [x(t) + /j, 3 ] is of the order of 1 un- 
less C < 1. Note that D c is scaled with N/MQ 2 , which 
is approximately equal to the slip distance. We wish to 
stress that, in granular layers, the critical slip distance is 
not scaled with the characteristic length of the microscopic 
geometry such as the surface roughness. This is the natu- 
ral consequence of the legth-free nature of the relaxation, 
represented in the form of equation (3). 

We also test a more plausible law for the dynamic friction, 
which is recently found in the DEM simulation. 

pss^Mo+s/"*, (14) 

where <f> ~ 0.3 and s is a numerical factor of the ord er of 
0.1 [Hatano, 2007]. Here we interpret I as V^m/Nd. The 
resulting dynamics is shown in Figure 4, where we adopt 
fio — 0.3, s = 0.2, and C = 1. The dynamics is very similar 
to that obtained in an experiment [Nasuno et al., 1998], in 
which the nondimensional number C is estimated to be of 
the order of 1. 

Despite the feasibility in reproducing an experimental re- 
sult on granular matter, we have to remark that the nondi- 
mensional parameter C = fir may be very small under 
a seismogenic condition. For example, if we assume that 
P — 100 MPa, d = 10 /jm, and H — 1 cm, using equation 
(6), the relaxation time r is of the order of 10~ 4 s. As the 
seismic slip takes place in seconds, C is of the order of 10~ 4 
so that the critical slip distance is negligible compared with 
the total slip distance. However, note that the framework of 
the length-free evolution law is not limited to the relaxation 
process of the velocity profile. It is straightforward to ex- 
tend the present evolution law to incorporate any processes 
that affect the frictional strength. One of the most illustrat- 
ing examples is a mechanochemical effect such as thermal 
decomposition of calcite, the rate constant of which is on 
the order of 1 sec or even much larger depending on the 
temperature [TTirono et al. , 2007]. A plausible modeling is 
in progress and the result will be presented elsewhere. 
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